#########################################################
# Application: EU legislature
# replication Rasmussen 2010
#########################################################
# set working directory
# setwd()

# load source function
source("replication_Rasmussen/cbq_binary_function.R")
# load library
library(readstata13)
# read in data
dat_r <- read.dta13("replication_Rasmussen/EUP388675-supplementary_material.dta") # downloaded from http://journals.sagepub.com/doi/suppl/10.1177/1465116510388675
dat_r$var15 <- as.numeric(as.character(dat_r$var15))
dat_r$var5 <- as.integer(dat_r$var5)
dat_r$var5[which(dat_r$var5==1)] <- 0
dat_r$var5[which(dat_r$var5==2)] <- 1
dat_r$var16 <- as.integer(dat_r$var16)
dat_r$var16[which(dat_r$var16==1)] <- 0
dat_r$var16[which(dat_r$var16==2)] <- 1
dat_r$var10 <- as.integer(dat_r$var10)
dat_r$var10 <- dat_r$var10 - 1 
dat_r$var15 <- dat_r$var15 - 1998
dat_r1 <- subset(dat_r,dat_r$var13==0)

# Covariates: 
# 
# Policy distance rapporteur-EP median
# Big party group
# Country coherence of key negotiators
# Closeness of working relationship (calendar year)
# Urgency of legislation: existing legislation expires; No existing legislation
# Workload (pending legislation during first reading)
# Issue areas consulted (no. of consultative committees)
# New act
# Directive
# Inter-institutional disagreement

# Focusing on risk of delegation
dat_r1 <- subset(dat_r1,select = c(var16,
                                   var6,
                                   var3, 
                                   var5,
                                   var10,
                                   var15,
                                   var1,
                                   var4,
                                   var11,
                                   var14,
                                   var7,
                                   var17)
                 )
dat_r1 <- dat_r1[complete.cases(dat_r1),]

# Replicate findings in Model III of Table 2
model1 <- glm(var16 ~ var6 +var3+ var5 + var10+ var15+ factor(var1)+ var4+ var11+ var14+ var7+ var17,
              data=dat_r1,
              family = binomial(link="logit")
             )
pred <- predict(model1,dat_r1,type="response")
pred_y <- ifelse(pred>0.5,1,0)
pred_ratio <- sum(pred_y == dat_r1$var16)/385 # 76.6%

###############################
# counterfactuals
dat_r2 <- with(dat_r1, data.frame(var6 = seq(0,2,length.out = 1000), 
                                  var3=seq(0,2,length.out = 1000), 
                                  var5=1,
                                  var10=median(var10), 
                                  var15=median(var15), 
                                  var1="Yes", 
                                  var4=median(var4,na.rm=T), 
                                  var11=median(var11, na.rm=T), 
                                  var14="Yes", 
                                  var7="Yes", 
                                  var17="Yes"
                                  )
               )

preds <- predict(model1,dat_r2,type="response", se.fit=TRUE)

predm <- preds$fit
predl <- preds$fit - (1.96*preds$se.fit) # lower bounds
predu <- preds$fit + (1.96*preds$se.fit) # upper bounds

